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ABSTRACT 

Initial steps in the application of a low-order 
pane! method computational fluid dynamics (CFD) 
code to the calculation of aircraft dynamic stability and 
control (S&C) derivatives are documented. Several 
capabilities, unique to CFD but not unique to this 
particular demonstration, are identified and 
demonstrated in this paper. These unique capabilities 
complement conventional S&C techniques and they 
include the ability to: 1) perform maneuvers without 
the flow/kinematic restrictions and support interference 
commonly associated with experimental S&C 
facilities, 2) easily simulate advanced S&C testing 
techniques. 3) compute exact S&C derivatives with 
uncertainty propagation bounds, and 4) alter the flow 
physics associated with a particular testing technique 
from those observed in a wind or water tunnel test in 
order to isolate effects. Also presented are discussions 
about some computational issues associated with the 
simulation of S&C tests and selected results from 
numerous surface grid resolution studies performed 
during the course of the study. 

INTRODUCTION 

Aircraft stability and control (S&C) 
derivatives quantify the aerodynamic forces and 
moments used to model aircraft dynamics. S&C 
derivatives are used to predict, for example, the 
longitudinal short period, lateral pure roll, lateral Dutch 
roll, spin behaviors, and handling qualities sensed by 
pilots for a given configuration. Historically, aircraft 
static and dynamic S&C derivatives have been 
determined through wind tunnel tests, water tunnel 
tests, or with empirical formulations 1 8 based on prior 
tests. However, most of the wind tunnels and all of the 
water tunnels used for S&C tests operate at low free 
stream velocities and Reynolds numbers. These 
limitations restrict the range of flight conditions that 
can be adequately simulated. 

Many S&C testing facilities provide unique 
dynamic capabilities that are not typically available in 
the wind tunnels used for aircraft performance tests. 
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Because of their dynamic capabilities, these S&C 
facilities may be expensive to build, operate, and/or 
maintain. The ran^e of motions that can be performed 
in such facilities may be limited by the wind tunnel 
walls or by the test rig’s kinematic and vibrational 
restrictions. In addition, many conventional dynamic 
test facilities do not provide the capability to determine 
damping and cross derivatives individually; instead, 
these facilities can only provide measurements of 
combined derivatives. The measurement of these 
quantities in pairs is not the preferred situation; it is 
simply a kinematic reality of the facilities available. A 
companion paper 10 details several computational 
techniques to separate such paired dynamic derivatives. 
Some facilities are not equipped with a “slip ring” 
capability and thus, cannot perform continuous rotary 
motions; instead they rely on oscillatory motions to 
provide brief periods of steady rotational motion. 

The difficulty of many experimental S&C 
facilities to adequately model the flight characteristics 
of today's aircraft over the entire flight envelope, the 
cost of developing better S&C testing facilities, and the 
recent advent of novel morphing aircraft 
configurations 1113 all contribute to the difficulty of 
predicting S&C derivatives. Prior research efforts 1 ' 29 
using CFD to predict S&C derivatives have made some 
progress toward bridging the capability gap, but these 
efforts are very time consuming 15 and have failed to 
make a significant impact in the day to day business of 
S&C prediction. Because both CFD and experimental 
studies 30 31 have their unique disadvantages, aerospace 
companies spend millions of dollars fixing S&C 
problems discovered during certification flight tests or 
production use. The current situation for S&C 
prediction is similar to the situation seen before CFD 
methods were widely used to accurately predict 
aerodynamic performance. 

CFD offers capabilities and techniques that 
complement experimental S&C testing techniques for 
obtaining S&C derivatives. With CFD, support and 

* Senior Research Scientist, Multidisciplinary 
Optimization Branch, M/S 159. Senior Member AIAA 
** Engineering Co-op Student. Multidisciplinary 
Optimization Branch, M/S 159, Student Member 
AIAA 


1 



AIAA 2004-0377 


wall interference effects can be either modeled or 
eliminated in order to assess their impact. Because the 
vehicle orientations and movements within the flow 
can be arbitrarily specified to obtain time-dependent 
solutions, CFD can easily simulate advanced S&C 
testing techniques (possibly including those for which 
no testing facility currently exists). By applying 
various differentiation techniques 12 * 16 ' 32 33 to CFD 
codes, researchers can obtain exact derivatives of the 
forces and moments throughout the time history of a 
maneuver. CFD can independently eliminate or refine 
various physics features within a flow to isolate 
particular effects. Finally, CFD can help explain the 
causes and types of separated flows affecting S&C 
prediction. 

In theory, applying CFD to S&C prediction is 
not only feasible, but straight forwardl; one would 
simply grid the vehicle of interest and perform the 
required calculations. In practice, many issues must be 
addressed before CFD can accurately predict S&C 
derivatives, including: 

• the cultural differences that exist between the 
CFD and S&C communities, evidenced by 
differences in language, notation, and 
expectations of the two groups 

• the significant increase in computer resources 
(execution time, memory, and disk space) 
required to address the aircraft S&C derivative 
needs 15 16 beyond the substantial resources 
already required for standard aircraft CFD 
performance calculations 

• the need for uncertainty quantification of S&C 
derivatives used in flight simulations and 
multidisciplinary design efforts and the further 
increase in computational resources 34 ' 38 
associated with computing S&C predictions 
with uncertainty included 

• the absence of a set of acceptable test cases 
that the CFD community must solve and the 
degree of accuracy required for CFD to 
become a credible alternative to wind and 
water tunnel testing 

To address these concerns, a team was 
recently formed at the National Aeronautics and Space 
Administration (NASA) Langley Research Center 
(LaRC). The Computation Methods for Stability and 
Control (COMSAC) team includes key researchers in 
the CFD and S&C communities devoted to developing 
and applying CFD methods for S&C needs. A recent 
workshop hosted by the COMSAC team devoted to 
this topic was well attended, and the topic gained broad 
support from the aerodynamic community. The 
COMSAC team is currently examining many research 
endeavors in this arena. 12 29 COMSAC' s mission is to: 
understand the causes and types of separated flows 
affecting S&C prediction: use this understanding to 


develop computational methods to predict the S&C 
derivatives of interest; and create a long-term plan for 
synergistic use of experimental facilities, testing 
techniques, and CFD methods to close the gap between 
CFD and experimental S&C activities. This paper is 
among the first to be published under this new effort. 

Results presented in this paper are relevant to 
furthering S&C applications of CFD. The results are 
computed by using a low-order panel method code on a 
simplified F-16XL fighter configuration. This work 
represents the first of many steps needed in the 
application of CFD to S&C problems. The 
configuration was chosen because it has a wide variety 
of measured data available. The panel code was 
chosen, among other reasons, for its quick run time, 
ease of use, extensive motion capabilities and its ability 
to represent three-dimensional bodies within time- 
dependent flows. 

In some ways, the use of a low-order panel 
method CFD code to compute the forces and moments 
of a highly swept configuration like the F-16XL at high 
angles of attack is a worst-case scenario. The vortical 
and viscous effects, which dominate the F-16XL 
flowfield, at high angles of attack, are not modeled in 
the code. The authors of this paper certainly do not 
advocate the uninformed use of a low-order panel 
method code for the computation of viscous, vortical, 
bursting-vortical, or separated flows. However, the 
code was used anyway because of the great potential to 
preliminary design that could result from a successful 
application. The focus of this paper is to demonstrate 
several computational techniques that advance the 
ability of CFD to compute S&C derivatives; this paper 
also raises a number of significant issues associated 
with S&C computations. The authors do advocate 
using the appropriate level of physics fidelity with 
respect to the problem posed, which can only be 
assessed by examining all available levels of fidelity. 
Producing quantifiable and justifiable error bounds on 
the accuracy of results is a critical area that must be 
addressed when using CFD to predict S&C derivatives. 
This paper demonstrates a grid resolution study and 
uncertainty propagation techniques that can be used as 
a model for providing error bounds in the future. Error 
bounds help identify the limitations of the current 
method; eventually, researchers will be able to 
compare the error bounds from all levels of fidelity and 
determine the required fidelity for any S&C 
application. 

TECHNICAL APPROACH 

This study involves the application of the 
automatic differentiation of FORTRAN ( ADIFOR) 
tool, version 2.0D, 40,41 to the Panel Method Ames 
Research Center (PMARC) code, 42 version 14.10. 
ADIFOR 2.0D was applied to PMARC by using a wide 
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variety of input variables as potential independent 
variables of differentiation. The PMARC body and 
wind axis forces and moments were used as the 
dependent variables of differentiation. Since 
uncertainty estimates of all S&C parameters (forces 
and moments and their derivatives) are normally 
requested by most control law designers, the 
computation of both first and second derivatives with 
respect to selected independent variables of 
differentiation are demonstrated; the latter are used to 
enable the propagation of known or assumed input 
variable uncertainties through the code to determine the 
uncertainty effects on the computed forces and 
moments. The relationship for uncertainty 
propagation 36 ' 38 as a function of normally distributed, 
random variables is 




i=l 


dx. 


,j = l*n 


m = number of random inputs 
n = number of output functions 
<7 i = standard deviation of 
the random inputs 


( 1 ) 


In this case, the output function F may be any of the 
computed body or wind axis forces and moments, or 
their first derivatives with respect to parameters such as 
the angle of attack, angle of sideslip, and the steady 
rotational rates. Original and ADIFOR-generated 
PMARC code versions were executed for an array of 
oscillatory and rotary aircraft motions to compare with 
data generated for a 10% scale F-16XL model tested 9 
in the NASA LaRC 14- by 22-Foot Subsonic Tunnel. 
Other data from wind tunnel (18% scale model) 43 and 
water tunnel (2.5% scale model) 44 ' 45 tests is also 
available, but was largely ignored for this study. The 
computational F-16XL geometry, as shown in Figure 1, 
was modified from that of the actual aircraft. The 
engine inlet was faired over with a smooth, bulbous 
surface and some detail near the aft end of the aircraft 
was simplified; this same geometry was used in the 
water tunnel tests. Limited comparisons with 
experimental data 9 ' 43 are presented in this paper. Other 
configurations, including a generic transport wing- 
body, a Boeing 747 wing-body-pylon-nacelle 
configuration, and an F/A-I8 A/B fighter were also 
executed during the course of the present studies; the 
results of these studies may be the subject of a follow- 
on publication. Only F-16XL data is shown in this 
paper. 


Computational Fluid Dynamics Code 
The demonstrations presented here employ the 
PMARC code, 42 a low-order panel method CFD code 
that was modified by applying AD to enable efficient 
computation of exact first and second force and 
moment derivatives with respect to a wide variety of 
code inputs. This panel method CFD code was chosen 
for its fast execution and moderate computational 
resource requirements, ease of use. ability to simulate a 
wide variety of steady and time-dependent rotary and 
oscillatory aircraft motions, and its amenability to 
processing by AD to compute derivatives that are exact 
to the formulation accuracy of the code solution. 

Second derivatives have been included to allow for 
input variable uncertainty propagation through the code 
with respect to the output forces and moments and their 
first derivatives. From the combination of the above 
mathematical techniques and the panel method choice, 
a demonstration of the computation of aircraft dynamic 
S&C derivatives, with propagated uncertainty bounds, 
could be rapidly performed. Sample results are 
presented for the F-16XL fighter configuration 
modeled with a coarse surface grid, shown in three- 
view in Figure 1. 

The PMARC, 42 version 14.10 potential flow 
solver is a FORTRAN 77 code that can compute 
surface pressures, forces, and moments of arbitrary 
shapes. The code assumes inviscid, irrotational, and 
incompressible flow. Boundary layer and compressible 
corrections are available but are not implemented in 
this study. PMARC can also compute solutions of 
unsteady, time-varying flow conditions for a wide 
variety of aircraft motions. The original PMARC code 
has been modified to allow for the input of a free 
stream Mach number, angle of attack, and angle of 
sideslip in degrees, rather than simply the input of free 
stream velocity components. The code was also 
modified to allow ADIFOR processing by replacing the 
use of scratch I/O files with common block access, and 
by inserting a few lines of code into the program to 
identify independent and dependent variables of 
differentiation. 

The input file to the program uses a set of grid 
points describing the shape of the geometry as a set of 
panels. In this study, both the right and left halves of a 
modified F-16XL configuration are modeled in the 
PMARC input geometry file. The gray-scale in Figure 
1 reflects different surface patches in the PMARC 
geometry input file, however these patches cannot be 
readily identified from the figure. The surface 
resolution, with 984 surface points and 566 surface 
panels, is moderately coarse by CFD standards, yet the 
vehicle is still clearly identifiable. Higher surface 
resolution grid files were available and used. Several 
such studies were performed and are discussed in the 
sections “Grid Resolution Studies" and “Results”. 
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PMARC allows the user to describe only half the 
aircraft and mirror the solution in the x-z plane. 
However, this technique does not capture the effects of 
a nonzero angle of sideslip, which are of interest in the 
S&C arena. The input fde also specifies the flight 
condition, certain algorithmic parameters, and the user- 
defined position of the reference point about which all 
moments are summed. The forces and moments are 
also nondimensionalized with a user-specified 
reference area, mean aerodynamic chord length, and 
wingspan. 

The PMARC code includes a wide variety of 
possible aircraft motions: 

• straight translation with imposed angles 
of attack ( a ) and sideslip ( /? ) 

• pure rotary motions about the three axes 
described as constant rotational rates (p, 
q, and r) 

• planar oscillatory motions in both 
translation and rotation about each of 
three axes, each described by an 
amplitude and frequency 

• coning motions (illustrated subsequently) 
used within advanced experimental S&C 
facilities 

In contrast to the limitations of some experimental 
facilities, coning motions can be simulated within 
PMARC with either continuous (i.e., a "slip ring") or 
limited roll motion capability (oscillatory test motion). 

In addition, the authors also use two other 
capabilities of the CFD code that are not generally 
available to S&C experimentalists. These capabilities 
are: the choice of either rigid or flexible wakes, and the 
ability to somewhat arbitrarily prescribe the time 
stepping characteristics of the solution. Although 
flexible wakes are preferred because they more 
accurately model the flow physics involved, the rigid 
wakes were more robust for some of the motions of 
interest and also lend themselves better to flow 
visualization. Wakes are currently attached to the wing 
trailing edges and the upper vertical tail trailing edge. 
The time stepping characteristics of a given solution 
can be defined to match the data sampling rate of a 
particular S&C facility, if desired. A caveat of this 
technique is that the size of the shed wake panels in the 
PMARC solution is proportional to the chosen time 
step, which alters the flow physics to some extent if 
sufficiently small time steps are not used. 

Automatic Differentiation 

Automatic differentiation 40 41 ' 46 (AD) is a 
technique for derivative computation. AD relies on the 
fact that every function, no matter how complicated, is 
executed on a computer as a (potentially very long) 
sequence of elementary arithmetic operations and 


functions evaluations (i.e., sine or cosine). By 
repeatedly applying the chain rule of differential 
calculus to the composition of those elementary 
operations, derivative information can be computed 
exactly and automatically. 

The two approaches for computing derivatives 
with AD are the forward mode and the reverse mode. 
The forward mode applies the chain rule of 
differentiation to propagate, equation by equation, 
derivatives of intermediate variables with respect to the 
input variables. In contrast, the reverse (adjoint) mode 
propagates, in reverse through the program, the 
derivatives of the output variables with respect to the 
input variables. The forward mode is belter suited to 
problems with fewer input variables than output 
variables, whereas the reverse mode is better suited to 
problems with fewer output variables than input 
variables. Many hybrids of the forward and reverse 
modes are possible, with complementary tradeoffs in 
required random access memory (RAM), disk space, 
and execution time. 

The forward method of AD is implemented in 
the ADIFOR, version 2.0D 40 41 tool. The reverse mode 
and some second derivative capability are also 
implemented in the ADIFOR, version 3.0. 4(> Both tools 
were developed jointly by the Center for Research on 
Parallel Computation at Rice University and the 
Mathematics and Computer Sciences Division at 
Argonne National Laboratory. Both tools have been 
used with the PMARC code, but only results using the 
ADIFOR 2.0D tool are shown in this paper. In general, 
to apply ADIFOR to a given FORTRAN 77 code, the 
user is only required to specify the independent and 
dependent variable names of the target differentiation 
and a differentiation starting point in the program. The 
AD tool then determines the variables that require 
derivative computations. It formulates the appropriate 
forward or reverse mode derivative expressions for 
these variables and generates new FORTRAN 77 code 
for the computation of both the original simulation and 
the associated derivatives. The ADIFOR-generated 
code is compiled (with special libraries) and executed 
like the original code. 

The forward mode of differentiation was used 
to compute the S&C derivatives because generally, 
only a few input, or independent, variables are 
identified for differentiation within a given code 
generation (compared with 12 output, or dependent, 
variables). The reverse mode of AD may be employed 
to calculate the S&C derivatives for morphing vehicles, 
as in References 12 and 13, where the potentially 
thousands of independent variables (perhaps, each 
surface grid point, along with the 30 or so flight 
condition and orientation variables) would greatly 
outnumber the 12 dependent variables. 


4 



AIAA 2004-0377 


Grid Resolution Studies 

Three assumptions were made at the start of 
this work: 1 ) a panel method CFD eode with a coarse 
surface grid could he used in conceptual or preliminary 
design whereas Euler / Navier-Stokes computations are 
likely to require so much execution time as to be 
incompatible with conceptual or preliminary design, 

2) the surface grid shown in Figure 1 was probably too 
coarse to be accurate but represented a good starting 
point, and 3) grid resolution could be increased at any 
time, once the computational processes were 
established, to smoothly increase the accuracy of the 
results. It remains to be seen if the first assumption is 
valid. However, as illustrated in Figure 2 of the 
Results section of the paper, the second assumption is 
valid; this surface grid inaccurately approximates even 
the most basic of S&C parameters, such as the lift 
variation with angle of attack, even at low angles of 
attack where a panel method formulation should be 
valid and vertical was not expected to affect the 
computation. This led the authors to, among other 
things, systematically examine the surface grid 
resolution requirements through a variety of studies, 
essentially testing the validity of the third assumption. 

A considerable amount of processing and 
analysis was involved in obtaining the surface grids to 
be used for the current surface grid resolution studies, 
which were obtained as subsets of volume grids. The 
F-16XL volume grid history is as follows: 

1 . The first grid was a 30-block, structured, 
half-configuration volume grid containing 
1,502,138 grid points (FINEST_VOLl); 
it was taken from previous Euler studies 
of References 39 and 43. Preliminary 
analysis using an earlier version of 
PMARC and an Euler/Navier-Stokes 
code was performed with this grid, but it 
provided unsuitable convergence for 
Navier-Stokes solutions in studies related 
to References 1 5 and 1 6 and was 
therefore abandoned. 

2. A new 36-block, structured, full 
configuration volume grid containing 
6,293,908 grid points (FINEST, VOL2) 
was developed for a modified F-16XL. 

An engine inlet fairing was added and tip 
missiles and launchers were removed. 
This configuration had improved grid 
resolution above the wing upper surface 
to enable users to capture votical flow 
over the upper wing surface with 
Euler/Navier-Stokes computations. 
Studies using this grid have not yet been 
performed. 


3. The Y -Z coordinate directions of the 
FlNEST_VOL2 grid were swapped to 
conform with PMARC requirements. 

4. Alternate grid points in each of the 1-J-K 
coordinate directions were removed to 
produce another volume grid 
(MED1UMJVOL, details presented in 
Table 1 ). 

5. Again, alternate grid points in each of the 
I-J-K coordinate directions were removed 
to produce yet another volume grid 
(COARSEJVOL). 

Surface grids consisting of 56 patches were extracted 
as subsets from each of the volume grids as follows: 

1 . A surface grid containing 1 1 3,856 grid 
points and 108,736 grid panels 
(F1NEST_SRF) was obtained from the 
FINEST, VOL2 grid. 

2. A reduced surface grid containing 
MEDIUM_SRF was obtained from the 
MEDIUM JVOL grid. 

3. A further reduced surface grid 
(COARSE_SRF) was obtained from the 
COARSEJVOL grid. 

4. Individual I-and J-grid lines from the 
COARSE_SRF grid were selected to be 
retained through two further sessions of 
coarsening (COARSEJSRF2 and 
COARSE_SRF3, respectively). 

All the above surface grids required some I-J 
reordering of grid lines to conform with the PMARC 
geometry convention for “correct side out” patches. 

The majority of runs performed during this study used 
the COARSE_SRF3, I-J surface grid, or variations of 
this grid (described subsequently). Details about all the 
volume and surface grids are presented in Table 1. 

Six types of surface grid resolution studies 
were conducted with the PMARC code. The studies 
used: 

1 . Externally generated grids for a 
symmetric, unswept, untapered, parabolic 
arc wing with an aspect ratio of two; 
alpha = 10 degrees (37 cases executed). 

2. Externally generated grids for 
configurations similar to the F-16XL 
wing plus fuselage carry-through 
(COARSE3_SRF) with parabolic arc 
wing sections; alpha = 10 degrees (45 
cases executed). 

3. Internally generated grids based on the 
original F-16XL configuration 
(COARSE3_SRF) input to PMARC, 
modified through input choices to 
PMARC; alpha = 10 degrees (34 cases 
executed). 
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4. Externally generated grids starting from 
higher resolution surface geometries 
(FINEST_SRF, MEDIUM_SRF, and 
COARSE_SRF): alpha = 10 degrees (3 
cases executed or attempted). 

5. Externally generated grids for 
configurations similar to the F-16XL 
wing plus part of the fuselage with 
parabolic arc wing sections; alpha = 30 
degrees (22 cases executed). 

6. Externally generated grids for 
configurations similar to the F-16XL 
wing plus part of the fuselage with 
parabolic arc wing sections; small 
amplitude forced pitch oscillation about 
alpha = 0 degrees (30 cases executed). 

The cases in group 4 failed to produce any useful 
results; cases in groups 5 and 6 produced results similar 
to the cases in group 2. Only the results from groups 1, 
2, and 3 are reported in this paper. The surface and 
volume grid history/usage is summarized in Table 1. 

RESULTS 

As mentioned previously. Figure 1 shows a 
three-view of the computational geometry of the 
F-16XL used for most of these calculations; the gray- 
scale indicate different surface patches in the PMARC 
geometry COARSE3_SRF input file. The surface 
resolution, consisting of 984 surface points and 566 
surface panels, is moderately coarse by CFD standards. 
Results are first presented for the static forces and 
moments, then for dynamic situations. Figure 2 shows 
a comparison of computed lift, pitching moment, and 
normal force coefficient data with measured data 9 43 
using this geometry. Circles in the figure are the lift 
force coefficient, triangles are the pitching moment 
coefficient, and squares are normal force coefficient; 
open symbols represent the computation and solid 
symbols are the data from Reference 9, as well as the 
solid diamonds representing the lift force coefficient 
from Reference 43. The computed data reflects the 
general character of the measured data, but the 
computed data differs significantly with both sets of 
measured data. It was assumed that most of the 
differences in the lift and normal force coefficient 
computed/measured data comparisons could be 
attributed to: 1 ) at high angles of attack, the inability of 
the PMARC code to model the vortical and viscous 
flows, 2) uncertainty propagation due to uncertainties 
in the angle of attack, and 3) the surface resolution 
used in the PMARC computations for all angles of 
attack. Minor differences between the tested geometry 
and the computational geometry may also contribute to 
the differences in the forces and moments. 

Vortex-generated lift is known to occur over 
the wing, 39 and the PMARC code cannot model this 


flow feature. However, several attempts were made to 
simulate the effects of vortical flow by attaching 
additional wakes to the upper wing surface near the 
leading edge and/or over the mid-chord of the wing. 
Two such attempts are shown in Figure 3; the figure 
includes the same lift force coefficient data as in Figure 
2, but with two additional curves representing the 
attempts to add additional wake to the wing upper 
surface, indicated with dashed lines. Indeed, the data 
comparisons could be improved at high angles of 
attack (30 to 50 degrees) by adding wakes separating 
from the mid-chord of the wing. However, such 
results required knowing either where the separation 
might occur or what value of the lift or normal force 
coefficient was to be obtained. The same technique 
significantly overpredicted the lift and normal force 
coefficient at 10 to 20 degrees angle of attack, 
indicated by the dashed curves extending off the top of 
the figure. This technique would need to be tailored for 
specific flow conditions. The cases should be repeated 
with both Euler and Navier-Stokes codes to truly assess 
the impact of vortex formation and bursting. 

The effect of uncertainty propagation through 
the CFD code was also examined: the same cases 
shown in Figure 2 were executed with an ADIFOR- 
generated Hessian (second derivative) code. For this 
case, derivatives of the forces and moments and their 
first derivatives with respect to angle of attack were 
differentiated again with respect to angle of attack. 
Uncertainties with respect to the angle of attack were 
computed and propagated through the code, 
proportional to variable gradients with respect to the 
uncertain variable, via Equation 1. The results from 
this study are shown as the error bars in Figures 3 and 
4 and in Table 2. Figure 3 shows the computed lift 
coefficient data of Figure 2, with added error bars 
showing the uncertainty in C/.due to a one-sigma 
uncertainty (for illustration purposes) in the angle of 
attack (AOA) for an assumed uncertainty in the 
measured angle of attack given by <J A0A = 0.1 degrees, 
determined from the measured data. Figure 4 shows 
the lift curve slope, C La , for the data in Figure 2, with 
added error bars showing uncertainty due to a one- 
sigma uncertainty in the angle of attack. Figures 3 and 
4 indicate that the propagated uncertainties due to 
uncertainties in angle of attack account for some, but 
only a small part, of the differences in the computed 
and measured lift coefficient. The data shown in Table 
2 illustrates a two-variable (angle of attack and pitch 
rate, Q) uncertainty assessment. The table shows the 
variable name, the uncertain variable name, the 
computed output variable sigma (from Eq. 1 ), and one-, 
two-, and three-sigma (representing 84.1. 97.7, and 
99.9% , respectively) low/high uncertainty bounds for 
each variable, as well as the mean value for each 
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variable. Table 2 illustrates a typical preliminary 
uncertainty analysis for the F-16XL configuration. The 

. . . _ dC, A dC. 

uncertainties in C, , — and — - (with a in 
L da dq 

degrees and q in degrees per second) are presented as 
calculated from first and second derivatives of C, with 
respect to a and q , assuming <7 = 0.1 for both a 

and q . The same capability is already available for all 
of the forces and moments in both the wind and body 
axis systems; results can be easily computed for any of 
me independeni variables of interest. 

The final element contributing to the 
discrepancy between measured and computed data in 
Figure 2 was the effect of surface resolution. A series 
of grid resolution studies, described previously, was 
performed. The results of these studies are shown in 
Figures 5 through 10. Figure 5 shows an envelope of 
computed lift coefficient data of various fidelities for a 
symmetric, unswept, untapered, parabolic arc wing 
with an aspect ratio of two, where alpha = 10 degrees 
(group I grid resolutions studies) at various surface 
resolutions, plotted as a function of the Log 10 (log 
base 10) of the number of surface panels. The 
envelope of solutions includes data from 37 cases 
(group 1 grid resolution studies) with the total surface 
(full span) resolution ranging from 8 to 4096 panels 
increasing in even multiples of 2, with distributions 
ranging from 2 to 256 panels chordwise, and 2 to 64 
panels spanwise, although not all possible 
combinations within that range were executed due to 
time and memory limitations on the available 
computers. Also shown in Figure 5 is an analytic 
computation of the lift coefficient taken from 
Reference 2. The envelope of computed lift coefficient 
data demonstrates grid independence at dense surface 
resolutions (thousands of surface panels); the data 
become both independent of the total number of 
surface panels and their distributions spanwise and 
chordwise, as well as agreeing well with the analytic 
solution. These results provide great evidence that the 
PMARC code can produce useful results for some 
geometries and flight conditions at sufficiently dense 
surface resolutions. 

Figure 6 illustrates the envelopes of data for 
the drag and the pitching moment coefficients. Similar 
to Figure 5, grid independence is observed in drag 
coefficient. The apparent lack of grid independence at 
dense surface resolutions in the pitching moment is 
simply an artifact of the magnified vertical scale 
compared with Figure 5. Nearly the same level of grid 
independence is achieved for both the lift force and 
pitching moment coefficient data. The convergence 
achieved in the drag coefficient is actually two orders 


of magnitude greater than either the lift force or 
pitehing moment coefficient. Figure 7 shows the same 
lift force coefficient data as in Figure 5, but plotted 
with curves of constant chordwise resolution (varying 
spanwise resolution) around the airfoil section. 
Comparing Figures 5 and 7, the minimum lift force 
boundary in Figure 5 is simply a composite of all the 
solutions with only four panels chordwise. In contrast, 
the maximum boundary in Figure 5 is a composite of 
the leftmost points on the curves in Figure 7 (those 
point with minimum spanwise resolution for a given 
chordwise resolution). A similar but inverse behavior 
was observed for the pitching moment. These results 
suggest, as expected, that four panels chordwise 
(around the airfoil) are too few; the resulting wing 
sections are modeled only as a composite of forward 
and rearward facing wedges, resulting in large errors in 
the predicted forces and moments. More panels 
chordwise appear to smoothly and significantly 
improve the fidelity of the PMARC solution. Put 
another way, the results only with only four panels 
chordwise never asymptotically approach the expected 
value with increased surface resolution for the cases 
studied. A discontinuous approach to the analytic 
solution is observed: the surface grid resolution must 
be "good enough” to begin an asymptotic approach to 
analytic value. A similar conclusion was reported in 
Reference 37. The F-16XL grid used for most of the 
studies has 10 panels chordwise around the airfoil 
section and 8 panels spanwise for each wing semispan; 
Based upon this result alone, one might expect the grid 
shown in Figure 1 to be sufficient to accurately 
approximate the measured results, which not the case. 

Figure 8 shows an envelope of lift coefficient 
data for wing shapes similar to the F-16XL planform 
(group 2 grid resolution studies, 45 cases). Behavior 
similar to that seen in Figure 5 is observed, but grid 
independence is not achieved to the same degree. The 
drag force and pitching moment data envelopes for 
these cases are similar to those of Figures 6. More 
panels, above a certain minimum threshold chordwise, 
improve the accuracy of the computed data. However, 
an important limiting effect in accuracy due to the 
more complex F-16XL wing geometry, compared to 
the swept tapered wing of Figures 5-7. is observed in 
Figure 8. 

Figure 9 shows a scatter plot (rather than an 
envelope of data) for the 34 cases of group 3 grid 
resolution studies. In this study, the total number of 
surface panels ranged from 504 to 2486, with the 
nominal grid containing 566 panels. The distribution 
of the panels varied from run to run, with surface 
resolution increasing either over the whole vehicle, 
only the wings, or select regions of the wings. This was 
a systematic study performed with the goal of matching 
the lift coefficient at 10 degrees angle of attack from 
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Reference 9. Many executions beyond those 34 
reported here were performed that resulted in true 
"outlier solutions” with unreasonably large lift 
coefficients (greater than 2.0); the same can also be 
said of the previous two grid resolution studies, 
although more "outlier solutions” occurred is the 
studies associated with Figure 9. Also shown in Figure 
9 are the values of the lift coefficient obtained from 
References 9 and 54. Most of the computed results in 
Figure 9 fall outside of the range of the two reference 
solutions, and no grid independence, nor any approach 
to the expected value, is discernable. Taken with 
previous results, this suggests that surface grid 
smoothness, in addition to the surface grid resolution 
and distribution, plays a big role in the accuracy of the 
PMARC results for lilt coefficient. In contrast, as 
shown in Figure 10, nearly any grid distribution or 
resolution may predict the expected pitching moment 
for this flow condition to acceptable accuracy. 

To summarize, several important effects have 
been identified in Figures 5-10 which strongly affect 
the accuracy of the PMARC results: 1) the total grid 
density, 2) the grid distribution spanwise and 
chordwise, 3) minimum thresholds of surface 
resolution may exist for obtaining accurate results, 4) 
configuration geometric complexity, and 5) grid 
smoothness. It is important to realize that once a grid 
has been obtained for a given configuration within a 
given study, the user may have little, if any, means to 
improve upon any of these aspects of the grid, which 
significantly affect the quality of the results. At this 
point, it is assumed that a grid of sufficient quality can 
be obtained and validated for use. The remainder of 
the figures illustrate in a qualitative sense what could 
be done with such a grid for S&C applications. 

Figure 1 1 shows the original vehicle 
geometry, with wake trailing from the top of the 
vertical tail, during a body axis roll maneuver. For all 
computations, other wakes (not shown in Figure 1 1) 
were also attached to the trailing edges of both wings. 
Figure 12 shows the vehicle, again with wake trailing 
from the top of the vertical tail, during a forced 
oscillation in roll maneuver. The reader should note 
the two reversals in rotary direction illustrated by the 
wake. In this case, the only useful “pure rotary” 
portion of the motion occurs between the two reversals 
after start-up transients have diminished; this is typical 
of how the maneuver might be performed in a wind 
tunnel test where the facility does not have “slip ring” 
capability. CFD offers the ability to perform the roll 
maneuver, without any support or wall interference 
effects, in either the pure rotary or in forced oscillation 
modes. 

Figures 13 and 14 illustrate two advanced 
experimental techniques used to determine dynamic 
S&C derivatives. Figure 13 shows a coning motion; 


the body rotates about the velocity vector, which is 
aligned with oncoming (low stream. Figure 14 
illustrates an oscillatory coning (or, inclined axis roll 
(IAR)) motion where the body rotates around an axis 
not aligned with the velocity vector. The motion is 
named for the oscillatory behavior of the angles of 
attack and sideslip observed during the motion. These 
types of motions can only be performed in a few wind 
and water tunnels globally. Figures 15 and 16 illustrate 
the CFD simulation oflhese two coning motions. In 
Figure 15, the original geometry is rotated to an initial 
pilch orientation angle, here 30.8 degrees, before the 
CFD solution is performed. A body axis roll maneuver 
is then performed with the rotated geometry and with 
the velocity vector aligned with the rotational vector, as 
above in Figure 1 1 . The motion could also be 
performed with an unrotated body and imposed angle 
of attack, but the body-axis forces differ between the 
two scenarios. In Figure 16, an angle of attack of 50 
degrees has been imposed on the body axis roll of the 
rotated geometry illustrated in Figure 15. The angle of 
attack in Figure 16, roll rates, and the forward speeds 
of Figures 11, 12, 15, and 16 were chosen simply to 
illustrate the maneuvers and do not reflect actual tunnel 
test conditions, although the body rotation angle of 
30.8 was chosen to match experimental data available 
for the F-16XL tested in the NASA LaRC 14- by 22- 
Foot Subsonic Tunnel. Figures 15 and 16 are included 
simply to illustrate the ease with which advanced 
testing techniques can be simulated through CFD. 
Obviously, the maneuvers should he performed with a 
grid of sufficient density, using guidance obtained from 
the above grid resolutions studies, and with a code that 
is able to capture all the relevant flow physics. 

Typical execution times for the F-16XL 
geometry in the original PMARC code are about 3.25 
minutes each on a R 12000 dual processor of a Silicon 
Graphics® Octane®. The forward mode ADIFOR- 
generated PMARC aerodynamic analysis, including the 
computation of angle of attack and sideslip derivatives, 
required approximately 6.29 minutes for the same 
single-processor execution scheme. The second 
derivative calculation, required for the uncertainty 
propagation analysis, required for the same two 
independent variables required about 14.10 minutes. 

All the calculations were performed in 64-bit 
arithmetic. The RAM requirements were about 1 30 
MB for the PMARC function and 2 1 3 MB for the 
function plus its derivatives with respect to angle of 
attack and angle of sideslip. The second derivative 
code required about 481 MB of storage. Late in the 
studies, the PMARC code was found to execute three 
to five times faster per case on a PC cluster, composed 
of a host and 32 worker nodes. Each worker node uses 
a 2193 Mhz central processing unit (CPU), the Linux 


8 



A1AA 2004-0377 


2.4.18 operating system, and a Lahey/Fujitsu 
Optimizing Fortran 95 compiler. 

CONCLUSIONS 

The Panel Method Ames Research Center 
(PM ARC) code, version 14.10, has been differentiated 
in several versions to provide exact first and second 
derivatives of the forces and moments with respect to 
wide variety of inputs. The original and differentiated 
codes have been executed to simulate a variety of 
motions typically used in experimental facilities to 
obtain dynamic stability and control (S&C) derivatives, 
initiai results, including input uncertainty propagation 
error bounds, using an F-16XL configuration as the test 
vehicle were shown. The results were obtained with a 
coarse surface grid resolution computational model 
using reasonable amounts of computer time and 
memory. The results generally reflect the trends of the 
measured F-16XL data, but the lift and normal force 
coefficient data were substantially different from the 
measured data. 

Additional wakes separating from the mid 
chord upper surface region of the wing could be used 
to improve the accuracy of the computed results, but 
this technique requires some knowledge about the 
expected flowfield characteristics. Propagated 
uncertainty bounds were also shown to account for a 
small part of the difference between the computed and 
measured data at low angles of attack. In addition, the 
results of numerous surface grid resolution studies 
were reported. The surface grid resolution studies 
show that a minimum of eight panels chordwise around 
the airfoil sections of a wing are necessary to obtain 
meaningful approximations to the expected forces and 
moments. The results are strongly influenced by the 
surface resolution in both the chordwise and spanwise 
directions, as well as the geometric complexity of the 
body and the smoothness of the grid distribution. For 
smooth distributions of surface panels, grid 
independence of results and agreement with known 
analytic solutions can be obtained with sufficiently 
dense surface resolutions. In the absence of smooth 
surface resolution, grid independence of the results 
cannot be achieved; in fact, the results in such cases do 
not appear to even converge to common values. 

Finally, the work reported here demonstrates 
some potential for using very coarse computational 
models to predict S&C characteristics of complex 
vehicles, especially when measured data is available to 
calibrate the computations. Such coarse computational 
models may be useful in the conceptual or preliminary 
design arena, but they represent only the first phase of 
applying CFD codes to aircraft S&C problems. The 
work also demonstrated several computational 
techniques that can be used to further the application of 
CFD to stability and control applications. 
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Table 1 Summary of F-16XL geometries, grids, and grid usages. 


Grid 

Description 

Grid Size 

Surface Resolution 


FINESTVOL1 

Half configuration 
volume grid w/tip 
missiles & launchers 

30 blocks 
1,502.138 grid 
points 

Not applicable 

Not used for the 
current studies 



36 blocks 
6.293.908 grid 
points 

Not applicable 

Future Euler and 
Navier-Stokes 
code use planned 

MEDIUM _VOL 

Subset of 

FINEST.VOL2 grid 

36 blocks 

837,924 grid points 

Not applicable 

Future Euler and 
Navier-Stokes 
code use planned 

COARSE. VOL 

SnsBBSI . 1 

36 blocks 

1 18,216 grid points 

Not applicable 

Future Euler and 
Navier-Stokes 
code use planned 

FINEST.SRF 

Surface grid 
extracted from 
FINEST.VOL2 
volume grid 

56 patches 
1 13,856 grid points 

108,736 surface panels 

Execution 
attempted, 
exceeded 
available RAM 

MEDIUM.SRF 

Surface grid 
extracted from 
MEDIUM.VOL 
volume grid 

56 patches 
29,772 grid points 

27,184 surface panels 

Execution 
attempted, 
exceeded 
available disk 
space 

COARSE.SRF1 

Surface grid 
extracted from 
CO ARSE_ V OL 
volume grid 

56 patches 
8,118 grid points 

6,796 surface panels 

Execution 
completed, 
unreasonably 
large forces and 
moments 
computed 

COARSE_SRF2 

Surface grid 
extracted from 
COARSE. VOL 
volume grid 

56 patches 
2,434 grid points 

1,744 surface panels 

Not used for the 
current studies 

COARSE.SRF3 

Surface grid 
extracted from 
CO ARSE. VOL 
volume grid 

56 patches 
984 grid points 

566 surface panels 

Used extensively 
for the current 
studies 


Table 2 Sample uncertainty analysis for F-16XL all values of a =0.1 . 



Variable Name 


dCL/dALPDEG 


dCL/dALPDEG 


dCL/dALPDEG 


CL 


dCL/dQ 


dCL/dQ 


dCL/dQ 


CL 


CL 


Uncertainty 

Variable 


ALPDEG 


Q 


ALPDEG+Q 


ALPDEG 


ALPDEG 


Q 


ALPDEG+Q 


ALPDEG 


TOTAL 


7TTTTT 


4.72E-05 


6.85E-05 




4.97E-05 


1 .23E-03 


1.23E-03 


2.41 E-03 


3.15E-03 


99 . 9 % 

Low 


2.02E-02 


2.02E-02 


2.02E-02 


3.95E-01 


2.39E-02 


2.04E-02 


2.04E-02 


3.94E-01 


3.91 E-01 


97 . 7 % 

Low 


2.03E-02 


2.03E-02 


2.02E-02 


3.97E-01 


2.40E-02 


2.16E-02 


2.16E-02 


3.96E-01 


3.95E-01 


84 . 1 % 

Low 


2.03E-02 


2.03E-02 


2.03E-02 


3.99E-01 


2.40E-02 


2.28E-02 


2.28E-02 


3.98E-01 


3.98E-01 


Mean 

Value 


2.04E-02 


2.04E-02 


EEEEBE1 


4.01 E-01 


2.41 E-02 


EHRBa 

M8M.a 


4.01 E-01 


2.04E-02 


2.04E-02 


2.04E-02 




2.41 E-02 


2.53E-02 


2.53E-02 


4.03E-01 


4.04E-01 


2.05E-02 


2.05E-02 


2.05E-02 


4.05E-01 


2.42E-02 


2.65E-02 


2.65E-02 


4.06E-01 


4.07E-01 


2.05E-02 


2.05E-02 


2.06E-02 


4.07E-01 


2.42E-02 


2.78E-02 


2.78E-02 


4.08E-01 


4.10E-01 
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Figure 1. Three-view of F-16XL computational model. 
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Figure 2. Computed static lift, moment, and 
normal force coefficient data and measured data. 



Figure 3. Computed static lift coefficient data, 
including uncertainty propagation error bars, and 
measured data. 



Figure 4. Computed static lift curve slope, C l/t , 
coefficient data, including uncertainty propagation 
error bars, and measured data. 



Figure 5. Envelope of computed lift coefficient data at 
various surface resolutions for AR = 2, symmetric, 
parabolic wing with analytic data. 



Figure 6. Envelopes of computed drag and pitching 
moment coefficient data at various surface resolutions, 
for AR = 2, symmetric, parabolic wing. 



Figure 7. Variation in computed moment coefficient 
with surface resolution for AR = 2, symmetric, 
parabolic wing and fixed chordwise resolutions. 



AIAA 



LoglO(Res) 

Figure 8. Envelope of computed lift coefficient data at 
various surface resolutions 
for configurations similar to F-16XL. 



LoglO(Res) 

Figure 9. Variation in computed lift coefficient with 
surface resolution 
for the F-16XL configuration. 



LoglO(Res) 

Figure 10. Variation in computed pitching moment 
coefficient with surface resolution 
for the F-16XL configurations. 
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Figure 1 1 . Sample body axis pure rotary roll 
maneuver, with trailing vertical tail wake. 



Figure 12. Sample body axis oscillatory roll maneuver, 
with trailing vertical tail wake. 




